__device__
float4 b2complex_new()
{
	return make_float4(0, 0, 0, 0);
}

__device__
float4 to_b2complex(float2 arg)
{
	return make_float4(arg.x/2, arg.y/2, 0, 0);
}

__device__
float2 to_complex(float4 arg)
{
	return make_float2(2*arg.y, 2*arg.w);
}

__device__
float4 b2complex_expi(float arg)
{
	return make_float4(cosf(arg)/2, cosf(arg)/2, sinf(arg)/2, sinf(arg)/2);
}

__device__
float4 b2complex_expj(float arg)
{
	return make_float4(cosf(arg)/2, cosf(arg)/2, -sinf(arg)/2, sinf(arg)/2);
}

__device__
float4 b2complex_add(float4 a, float4 b)
{
	return make_float4(a.x + b.x, a.y + b.y, a.z + b.z, a.w + b.w);
}

__device__
float4 b2complex_addu0u2(float4 a, float4 b)
{
	return make_float4(a.x + b.x, 0, a.z + b.z, 0);
}

__device__
float4 b2complex_addu1u3(float4 a, float4 b)
{
	return make_float4(0, a.y + b.y, 0, a.w + b.w);
}

__device__
float4 b2complex_sub(float4 a, float4 b){
	return make_float4(a.x - b.x, a.y - b.y, a.z - b.z, a.w - b.w);
}

__device__
float4 b2complex_subu0u2(float4 a, float4 b){
	return make_float4(a.x - b.x, 0, a.z - b.z, 0);
}

__device__
float4 b2complex_subu1u3(float4 a, float4 b){
	return make_float4(0, a.y - b.y, 0, a.w - b.w);
}

__device__
float4 b2complex_multyu0u2(float4 a, float4 b){
	return make_float4(	2*a.x*b.x - 2*a.z*b.z,
						0,
						2*a.x*b.z - 2*a.z*b.x,
						0);
}

__device__
float4 b2complex_multyu1u3(float4 a, float4 b){
	return make_float4(	0,
						2*a.y*b.y - 2*a.w*b.w,
						0,
						2*a.y*b.w + 2*a.w*b.y);
}

__device__
int mirror_index(int srcIndex, int matrixWidth)
{
	int matrixIndex = srcIndex/matrixWidth * matrixWidth;
	return (matrixIndex == srcIndex) ? srcIndex : matrixWidth - srcIndex + 2*matrixIndex;
}

__device__
void mirror(float4* src, int size, int down, int right, int dmirr, int rmirr)
{
	int downright = size*down + right;
	if (rmirr != right){
		int downrmirr = down*size + rmirr;
		src[downrmirr].x = src[downright].x;
		src[downrmirr].y = src[downright].y;
		src[downrmirr].z = -src[downright].z;
		src[downrmirr].w = -src[downright].w;
	}
	if (dmirr != down){
		int dmirrright = dmirr*size + right;
		src[dmirrright].x = src[downright].x;
		src[dmirrright].y = -src[downright].y;
		src[dmirrright].z = src[downright].z;
		src[dmirrright].w = -src[downright].w;
	}
	if (rmirr != right && dmirr != down){
		int dmirrrmirr = dmirr*size + rmirr;
		src[dmirrrmirr].x = src[downright].x;
		src[dmirrrmirr].y = -src[downright].y;
		src[dmirrrmirr].z = -src[downright].z;
		src[dmirrrmirr].w = src[downright].w;
	}
}


__device__ __constant__ float _PI = 3.141592653;

extern "C"
__global__ void fftu0u2(int size, float4* a, float4* wilist, float4* wjlist, int iteration){
	int idx = blockIdx.x * gridDim.x + threadIdx.x;
	int idy = blockIdx.y * gridDim.y + threadIdx.y;
	if (idx < size && idy < size){
		int step = 1 << iteration;
		int step2 = step*2;
		bool flag = (idy/step) % 2 + (idx/step) % 2 > 0;
		if (!flag){
			float4 w1 = wilist[(idy % step2)*(size/step2)];
			float4 w2 = wjlist[(idx % step2)*(size/step2)];
			
			int a00i = idy*size + idx;
			int a10i = (idy+step)*size + idx;
			int a01i = idy*size + (idx+step);
			int a11i = (idy+step)*size + (idx+step);
			float4 a00 = a[a00i];
			float4 a10 = a[a10i];
			float4 a01 = a[a01i];
			float4 a11 = a[a11i];
			
			float4 newa00 = b2complex_new();
			float4 newa01 = b2complex_new();
			float4 newa10 = b2complex_new();
			float4 newa11 = b2complex_new();
			
			newa00 = b2complex_addu0u2(newa00, a00);
			newa10 = b2complex_addu0u2(newa10, a00);
			newa01 = b2complex_addu0u2(newa01, a00);
			newa11 = b2complex_addu0u2(newa11, a00);
			
			float4 t = b2complex_multyu0u2(w1, a10);
			newa00 = b2complex_addu0u2(newa00, t);
			newa10 = b2complex_subu0u2(newa10, t);
			newa01 = b2complex_addu0u2(newa01, t);
			newa11 = b2complex_subu0u2(newa11, t);
			
			t = b2complex_multyu0u2(w2, a01);
			newa00 = b2complex_addu0u2(newa00, t);
			newa10 = b2complex_addu0u2(newa10, t);
			newa01 = b2complex_subu0u2(newa01, t);
			newa11 = b2complex_subu0u2(newa11, t);
			
			t = b2complex_multyu0u2(w2, a11);
			t = b2complex_multyu0u2(w1, t);
			newa00 = b2complex_addu0u2(newa00, t);
			newa10 = b2complex_subu0u2(newa10, t);
			newa01 = b2complex_subu0u2(newa01, t);
			newa11 = b2complex_addu0u2(newa11, t);
			
			
			a[a00i] = b2complex_add(a00, newa00);
			a[a10i] = b2complex_add(a10, newa10);
			a[a01i] = b2complex_add(a01, newa01);
			a[a11i] = b2complex_add(a11, newa11);
		}
	}
	__syncthreads();
}

extern "C"
__global__ void fft(int size, float4* a, float4* wilist, float4* wjlist, int iteration){
	int idx = blockIdx.x * gridDim.x + threadIdx.x;
	int idy = blockIdx.y * gridDim.y + threadIdx.y;
	if (idx < size && idy < size){
		int step = 1 << iteration;
		int step2 = step*2;
		bool flag = (idy/step) % 2 + (idx/step) % 2 > 0;
		if (!flag){
			float4 w1 = wilist[(idy % step2)*(size/step2)];
			float4 w2 = wjlist[(idx % step2)*(size/step2)];
			
			int a00i = idy*size + idx;
			int a10i = (idy+step)*size + idx;
			int a01i = idy*size + (idx+step);
			int a11i = (idy+step)*size + (idx+step);
			float4 a00 = a[a00i];
			float4 a10 = a[a10i];
			float4 a01 = a[a01i];
			float4 a11 = a[a11i];
			
			float4 newa00 = b2complex_new();
			float4 newa01 = b2complex_new();
			float4 newa10 = b2complex_new();
			float4 newa11 = b2complex_new();
			
			newa00 = b2complex_addu1u3(newa00, a00);
			newa10 = b2complex_addu1u3(newa10, a00);
			newa01 = b2complex_addu1u3(newa01, a00);
			newa11 = b2complex_addu1u3(newa11, a00);
			
			float4 t = b2complex_multyu1u3(w1, a10);
			newa00 = b2complex_addu1u3(newa00, t);
			newa10 = b2complex_subu1u3(newa10, t);
			newa01 = b2complex_addu1u3(newa01, t);
			newa11 = b2complex_subu1u3(newa11, t);
			
			t = b2complex_multyu1u3(w2, a01);
			newa00 = b2complex_addu1u3(newa00, t);
			newa10 = b2complex_addu1u3(newa10, t);
			newa01 = b2complex_subu1u3(newa01, t);
			newa11 = b2complex_subu1u3(newa11, t);
			
			t = b2complex_multyu1u3(w2, a11);
			t = b2complex_multyu1u3(w1, t);
			newa00 = b2complex_addu1u3(newa00, t);
			newa10 = b2complex_subu1u3(newa10, t);
			newa01 = b2complex_subu1u3(newa01, t);
			newa11 = b2complex_addu1u3(newa11, t);
			
			a[a00i] = newa00;
			a[a10i] = newa10;
			a[a01i] = newa01;
			a[a11i] = newa11;
		}
	}
	__syncthreads();
}

extern "C"
__global__ void fftMirror(int size, float4* a, float4* wilist, float4* wjlist, int iteration){
	int idx = blockIdx.x * gridDim.x + threadIdx.x;
	int idy = blockIdx.y * gridDim.y + threadIdx.y;
	/*if (idx < size && idy < size){
		int step = 1 << iteration;
		int matrixWidth = step << 1;
		int part = matrixWidth/4;
		int dmatrInd = idy/matrixWidth * matrixWidth;
		int rmatrInd = idx/matrixWidth * matrixWidth;
		int tdown = idy - dmatrInd;
		int tright = idx - rmatrInd;
		bool flag = (tdown <= part && tright <= part);
		if (flag){
			float4 w1 = wilist[(idy % matrixWidth)*(size/matrixWidth)];
			float4 w2 = wjlist[(idx % matrixWidth)*(size/matrixWidth)];
			
			int a00i = idy*size + idx;
			int a10i = (idy+step)*size + idx;
			int a01i = idy*size + (idx+step);
			int a11i = (idy+step)*size + (idx+step);
			float4 a00 = a[a00i];
			float4 a10 = a[a10i];
			float4 a01 = a[a01i];
			float4 a11 = a[a11i];
			
			float4 newa00 = ccomplex_new();
			float4 newa01 = ccomplex_new();
			float4 newa10 = ccomplex_new();
			float4 newa11 = ccomplex_new();
			
			newa00 = ccomplex_add(newa00, a00);
			newa10 = ccomplex_add(newa10, a00);
			newa01 = ccomplex_add(newa01, a00);
			newa11 = ccomplex_add(newa11, a00);
			
			float4 t = ccomplex_multy(w1, a10);
			newa00 = ccomplex_add(newa00, t);
			newa10 = ccomplex_sub(newa10, t);
			newa01 = ccomplex_add(newa01, t);
			newa11 = ccomplex_sub(newa11, t);
			
			t = ccomplex_multy(w2, a01);
			newa00 = ccomplex_add(newa00, t);
			newa10 = ccomplex_add(newa10, t);
			newa01 = ccomplex_sub(newa01, t);
			newa11 = ccomplex_sub(newa11, t);
			
			t = ccomplex_multy(w2, a11);
			t = ccomplex_multy(w1, t);
			newa00 = ccomplex_add(newa00, t);
			newa10 = ccomplex_sub(newa10, t);
			newa01 = ccomplex_sub(newa01, t);
			newa11 = ccomplex_add(newa11, t);
			
			a[a00i] = newa00;
			a[a10i] = newa10;
			a[a01i] = newa01;
			a[a11i] = newa11;
			
			if (iteration >= 2){
				int rmirr = mirror_index(idx, matrixWidth);
				int dmirr = mirror_index(idy, matrixWidth);
							
				mirror(a, size, idy, idx, dmirr, rmirr);
							
				rmirr = mirror_index(idx + step, matrixWidth);
				mirror(a, size, idy, idx + step, dmirr, rmirr);
							
				dmirr = mirror_index(idy + step, matrixWidth);
				rmirr = mirror_index(idx, matrixWidth);
				mirror(a, size, idy + step, idx, dmirr, rmirr);
							
				rmirr = mirror_index(idx + step, matrixWidth);
				mirror(a, size, idy + step, idx + step, dmirr, rmirr);
			}
		}
	}*/
	__syncthreads();
}

extern "C"
__global__ void wlist(int size, float4* wi, float4* wj){
	int blockSize = blockDim.x * blockDim.y;
	int block = blockIdx.y * gridDim.x + blockIdx.x;
	int idx = (blockSize * block) + threadIdx.y * blockDim.x + threadIdx.x;
	if (idx < size){
		wi[idx] = b2complex_expi(-2*_PI*idx/size);
		wj[idx] = b2complex_expj(-2*_PI*idx/size);
	}
}

extern "C"
__global__ void inverseRows(int size, float4* src){
	int blockSize = blockDim.x * blockDim.y;
	int block = blockIdx.y * gridDim.x + blockIdx.x;
	int idx = (blockSize * block) + threadIdx.y * blockDim.x + threadIdx.x;
	int col = idx;
	if (idx < size){
		int icol = 0;
		int jcol = 0;
		int dl = size/2;
		int st = size - 1;
		int j = 0;
		for (int i = 0; i < st; i++){
			if (i < j){
				icol = i*size + col;
				jcol = j*size + col;
				float4 s0 = src[icol];
				src[icol] = src[jcol];
				src[jcol] = s0;
			}
			int k = dl;
			while (k <= j){
				j = j - k;
				k = k >> 1;
			}
			j = j + k;
		}
	}
}

extern "C"
__global__ void inverseCols(int size, float4* src){
	int blockSize = blockDim.x * blockDim.y;
	int block = blockIdx.y * gridDim.x + blockIdx.x;
	int idx = (blockSize * block) + threadIdx.y * blockDim.x + threadIdx.x;
	int row = idx;
	if (idx < size){
		int rowi = 0;
		int rowj = 0;
		int dl = size/2;
		int st = size - 1;
		int j = 0;
		for (int i = 0; i < st; i++){
			if (i < j){
				rowi = row*size + i;
				rowj = row*size + j;
				float4 s0 = src[rowi];
				src[rowi] = src[rowj];
				src[rowj] = s0;
			}
			int k = dl;
			while (k <= j){
				j = j - k;
				k = k >> 1;
			}
			j = j + k;
		}
	}
}

extern "C"
__global__ void inverse(int size, int power, float4* src){
	int idx = blockIdx.x * gridDim.x + threadIdx.x;
	int idy = blockIdx.y * gridDim.y + threadIdx.y;
	if (idx < size && idy < size){
		
		int indexy = idy;
		int indexx = idx;
		int newrow = 0;
		int newcol = 0;
		for (int l = power; l >= 0; l--){
			int k = 1 << l;
			if (k <= indexy){
				indexy -= k;
				int delta = 1 << (power-l);
				newrow += delta;
			}
			if (k <= indexx){
				indexx -= k;
				int delta = 1 << (power-l);
				newcol += delta;
			}
		}
		int rowcol = idy*size + idx;
		int newrownewcol = newrow*size + newcol;
		if (idy > newrow){
			float4 t = src[rowcol];
			src[rowcol] = src[newrownewcol];
			src[newrownewcol] = t;
		} else if (idy == newrow && idx > newcol){
			float4 t = src[rowcol];
			src[rowcol] = src[newrownewcol];
			src[newrownewcol] = t;
		}
	}
}

extern "C"
__global__ void convertToCComplex(int size, float2* src, float4* out){
	int idx = blockIdx.x * gridDim.x + threadIdx.x;
	int idy = blockIdx.y * gridDim.y + threadIdx.y;
	if (idx < size && idy < size){
		int uid = idy*size + idx;
		out[uid] = to_b2complex(src[uid]);
	}
}

extern "C"
__global__ void convertToComplex(int size, float2* src, float4* out){
	int idx = blockIdx.x * gridDim.x + threadIdx.x;
	int idy = blockIdx.y * gridDim.y + threadIdx.y;
	if (idx < size && idy < size){
		int uid = idy*size + idx;
		src[uid] = to_complex(out[uid]);
	}
}
